library(tidyverse)
library(estimatr)
library(marginaleffects)
library(ggsci)
library(Cairo)
library(patchwork)

#csvデータを読み込む
data <- read_csv("datafile.csv")

#実験群
data <- data %>% 
  mutate(exp = 
           case_when(!is.na(con) ~ "統制群",
                     !is.na(populism) ~ "ポピュリズム群",
                     !is.na(conspiracy) ~ "陰謀論群"),
         exp = factor(exp, levels = c("統制群","ポピュリズム群","陰謀論群")),
         out1 = 6 - Q13.4,
         sat = 
           case_when(con == 3 ~ "ok",
                     populism == 3 ~ "ok",
                     conspiracy == 3 ~ "ok",
                     TRUE ~ "ng") )

#不良回答者を除く
data <- data %>% 
  filter(sat == "ok")


# 加重回帰（単純なウェイト付き最小二乗）
TukeyHSD(aov(data$out1 ~ data$exp))

#図1の分析
sum <- data %>% 
  group_by(exp) %>% 
  summarise(mean = mean(out1, na.rm = T),
            sd = sd(out1, na.rm = T),
            n = n(),
            se = sd / sqrt(n),
            ci = se * 1.96) %>% 
  mutate(exp_label = paste0(exp, "\n (n = ", n, ")")) %>% 
  mutate(exp_label = factor(exp_label, levels = c("統制群\n (n = 638)","ポピュリズム群\n (n = 649)","陰謀論群\n (n = 623)")))

table(sum$exp_label)

#図1の可視化
sig <- data.frame(x = c(1, 1, 2.1),
                  xend = c(2, 3, 3),
                  y = c(3.5, 3.8, 3.5),
                  yend = c(3.6, 3.9, 3.6),
                  label = c("p = .42", "p = .09", "p = .66"),
                  label.x = c(1.5, 2, 2.5),
                  label.y = c(3.7, 4, 3.7))

fig <- sum %>%
  mutate(xmin = as.numeric(factor(exp_label)) - 0.3,
         xmax = as.numeric(factor(exp_label)) + 0.3,
         ymin = 1,
         ymax = mean) %>%
  ggplot() +
  geom_rect(aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), fill = "gray50") +
  geom_errorbar(aes(x = as.numeric(factor(exp_label)), ymin = mean - ci, ymax = mean + ci),
                width = 0.2, size = 0.8) +
  geom_label(aes(x = as.numeric(factor(exp_label)), y = mean, label = sprintf("%.2f", mean)),
             vjust = 5, size = 5) +
  geom_text(aes(x = label.x, y = label.y, label = label), data = sig) +
  geom_segment(aes(x = x, xend = x, y = y , yend = yend), data = sig) +
  geom_segment(aes(x = x, xend = xend, y = yend, yend = yend), data = sig) +
  geom_segment(aes(x = xend, xend = xend, y = y, yend = yend), data = sig) +
  labs(x = "実験群", y = "平均値") +
  scale_x_continuous(breaks = 1:3, labels = sum$exp_label) +
  scale_y_continuous(breaks = seq(1, 4, 0.5)) +
  coord_cartesian(ylim = c(1, 4)) +
  theme_bw(base_size = 14, base_family = "HiraginoSans-W3") +
  theme(axis.title.y = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(margin = margin(t = 5, b = -5), size = 12)) +
  guides(size = "none")

fig

quartz(file = "単純統計.pdf", 
       type = "pdf", family = "HiraginoSans-W3",
       width = 9, height = 4)
print(fig)
dev.off()

###########################################################################################
#陰謀論的信念
data2 <- data %>% 
  drop_na(con1:con5)

#主成分分析
pca_result <- prcomp(data2 %>% 
                       dplyr::select(con1:con5),
                     center = TRUE, scale. = TRUE)

pca_loadings <- as.data.frame(pca_result$rotation) %>%
  rownames_to_column(var = "variable") 

#主成分分析の可視化
pca_fig <- pca_loadings %>% 
  mutate("第1主成分" = PC1 * -1,
         "第2主成分" = PC2* -1) %>% #主成分得点を逆転
  pivot_longer(- variable) %>% 
  filter(variable != "ID") %>% 
  mutate(variable = factor(variable, levels = c("con1","con2","con3","con4","con5"))) %>% 
  filter(name == "第1主成分" | name == "第2主成分")


fig0 <- pca_fig %>% 
  mutate(variable = fct_rev(variable))%>%
  ggplot() +
  geom_segment(aes(y = variable, yend = variable,
                   x = 0, xend = value)) +
  geom_point(aes(y = variable, x = value), color = "black", size = 3) +
  labs(x = "主成分負荷量", y = "陰謀論的信念に関する変数") +
  geom_vline(xintercept = 0, linetype = "solid", color = "grey") +
  theme_bw(base_size = 12, base_family = "HiraginoSans-W3") +
  theme(strip.text = element_text(size = 12),
        axis.ticks.y = element_blank(),
        axis.text.y = element_text(size = 10)) + 
  facet_wrap(~ name) + 
  geom_vline(xintercept = 0.35, linetype = "dashed", color = "gray40") + 
 # scale_x_continuous(breaks = seq(-.2, 1, by = .2),limits=c(-.2,1)) + 
  geom_text(aes(x = value, y = variable, label=sprintf("%.2f", value)),
            size = 3, nudge_y = .3)

fig0


###############################################################################
#主成分得点とマージ
data2 <- data2 %>%
  mutate(pc1 = pca_result$x[, 1] * -1)

#回帰分析
fit <- data2 %>% 
  lm_robust(out1 ~ exp * pc1, 
            data = .,
            se_type = "HC2")

#限界効果
newdata <- datagrid(model = fit,
                    exp = c("統制群", "ポピュリズム群", "陰謀論群"),
                    pc1 = seq(-4.94, 4.43, length.out = 100))
#95%信頼区間
box95 <- slopes(fit,
               variables = "exp",
               newdata = newdata,
               slope = "dydx",
               conf_level = 0.95) %>% 
  mutate(sig = if_else(conf.low * conf.high > 0, "sig", "ns"),
         team = "l95")

#90%信頼区間
box90 <- slopes(fit,
               variables = "exp",
               newdata = newdata,
               slope = "dydx",
               conf_level = 0.90) %>% 
  rename(cl = conf.low,
         ch = conf.high) %>% 
  dplyr::select(cl, ch)

#マージ
box1 <- cbind(box95,box90) %>% 
  mutate(contrast = 
           case_when(contrast == "ポピュリズム群 - 統制群" ~ "ポピュリズム群",
                     contrast == "陰謀論群 - 統制群" ~ "陰謀論群"))


#図2左      
fig1 <- box1 %>% 
  ggplot(aes(y = estimate, x = pc1)) +
  geom_line() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.25) +
  geom_ribbon(aes(ymin = cl, ymax = ch), alpha = 0.5) +
  facet_wrap(~ contrast) + 
  labs(x = "陰謀論的信念", y = "陰謀論的信念が実験処置の効果に与える影響",
       subtitle = "限界効果") + 
  theme_bw(base_size = 14, base_family = "HiraginoSans-W3") + 
  theme(axis.title.y = element_text(size = 12),
        axis.text.y = element_text(size = 10),
        legend.position = "bottom",
        axis.title.x = element_text(margin = margin(t = 5, b = -5), size = 12)) + 
  scale_y_continuous(breaks = seq(-1,1, by=.25),limits=c(-1, 1)) + 
  geom_hline(yintercept =  0, color = "black",  linetype="dashed")

fig1


###################################################################################

# 予測値
pred_fig <- predictions(fit,
                       newdata = datagrid(
                         exp = c("統制群","ポピュリズム群","陰謀論群"),
                         pc1 = c(-4.9,0,4.4)))

#図2右可視化
pred_fig <- pred_fig %>% 
  mutate(陰謀論的信念 = 
           case_when(pc1 == -4.9 ~ "低群",
                     pc1 == 0 ~ "中群",
                     pc1 == 4.4 ~ "高群"),
         陰謀論的信念 = factor(陰謀論的信念, levels = c("低群","中群","高群")))

fig2 <- pred_fig %>% 
  ggplot(aes(x = exp, y = estimate, shape = 陰謀論的信念)) +
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high),
                  position = position_dodge(.7), size = .8) +
  labs(x = "実験群", y = "予測値",shape = "陰謀論的信念",
       subtitle = "予測値") + 
  theme_bw(base_size = 14, base_family = "HiraginoSans-W3") + 
  theme(axis.title.y = element_text(size = 12),
        axis.text.y = element_text(size = 10),
        legend.position = "bottom",
        axis.title.x = element_text(margin = margin(t = 5, b = -5), size = 12)) + 
  geom_text(aes(x = exp, y = estimate, label=sprintf("%.2f", estimate), group = 陰謀論的信念),
            position = position_dodge(width = 0.7),
            hjust = 1.5)

fig2

#####################################################################
#限界効果の図と予測値の図をマージ
fig_23 <- fig1 | fig2 + 
  plot_annotation(tag_levels = 'A') 

fig_23

quartz(file = "限界効果＋予測値.pdf", 
       type = "pdf", family = "HiraginoSans-W3",
       width = 12, height = 5)
print(fig_23)
dev.off()

